// ZakharovFunction.cpp: implementation of the ZakharovFunction class.
// when run on 20 dimensions, the output in zakharov1.out is interesting
//////////////////////////////////////////////////////////////////////
#include <iomanip>
#include <fstream>
#include <string>

#include "zakharov_function.h"
#include "DataHolder.h"



//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

ZakharovFunction::ZakharovFunction(int i_NInput)
{
  NInput = i_NInput;

  // Create 2 arrays
  newA(mArray, NInput); 
  newA(mDerivative, NInput);

  Result.Size(NInput);
}


ZakharovFunction::~ZakharovFunction()
{
  deleteA(mArray);
  deleteA(mDerivative);
}

// takes a point and if it is beyong boundary move the point
// into the boundary box
void ZakharovFunction::checkBoundary(vector<double> &Input)
{
  for ( unsigned int i = 0; i < Input.size(); i++)
  {
    // reflective boundary
    /*while( (Input[i] >= 10.0) || (Input[i] <= -5.0) )
    {
      if ( Input[i] >= 10.0 )
      {
        Input[i] = 10.0 - fabs(Input[i]-10.0);
      }
      else if ( Input[i] <= -5.0 )
      {
        Input[i] = -5.0 + fabs(-5.0-Input[i]);
      }
    }*/

    // periodic boundary
    if ( Input[i] >= 10.0 )
    {
      Input[i] = -5.0 + fmod(Input[i],10.0);
    }
    else if ( Input[i] <= -5.0 )
    {
      Input[i] = 10.0 + fmod(Input[i],5.0);
    }

    // hard boundaries
    /*if ( Input[i] >= 10.0 )
    {
      Input[i] = 10.0;
    }
    else if ( Input[i] <= -5.0 )
    {
      Input[i] = -5.0;
    }*/
  }
}

// takes a point and if it is beyong boundary move the point
// into the boundary box
void ZakharovFunction::checkBoundary(double *Input)
{
  for (int i = 1; i <= NInput; i++)
  {
    // periodic boundary
    if ( Input[i] >= 10.0 )
    {
      Input[i] = -5.0 + fmod(Input[i],10.0);
    }
    else if ( Input[i] <= -5.0 )
    {
      Input[i] = 10.0 + fmod(Input[i],5.0);
    }
  }
}

//////////////////////////////////////////////////////////////////////
// functions
//////////////////////////////////////////////////////////////////////

double ZakharovFunction::Value(double *Input)
{
  // first check that point is not beyong any boundaries
  checkBoundary(Input);

  // Initialise vars answer and danswer (both local)
  double answer = 0.0;
  double danswer = 0.0;

  // Constructing the equation
  for (int m = 1; m <=NInput; m++)
  {
    answer += (Input[m]*Input[m]);
  }

  // Makes danswer equal to itself multiplied
  for(int n=1; n<=NInput; n++)
  {
    danswer += ((0.5*n)*Input[n]);
  }

  answer += (danswer*danswer);
  answer += (danswer*danswer*danswer*danswer);

  Result.Input(answer, Input);

  nCall++;

  return (answer); 
}

// Calculates/builds the function value
double ZakharovFunction::Value(vector<double> &Input)
{
  // first check that point is not beyong any boundaries
  checkBoundary(Input);

  // Initialise vars answer and danswer (both local)
  double answer = 0.0;
  double danswer = 0.0;

  // Constructing the equation
  for (int m = 0; m < NInput; m++)
  {
    answer += (Input[m]*Input[m]);
  }

  for(int n=0; n<NInput; n++)
  {
    danswer += ((0.5*(n+1))*Input[n]);
  }

  answer += (danswer*danswer);
  answer += (danswer*danswer*danswer*danswer);

  Result.Input(answer, Input);

  nCall++;

  return (answer); 
}

double ZakharovFunction::NormValue(double *Input)
{
  // first check that point is not beyong any boundaries
  checkBoundary(Input);

  // Initialise vars answer and danswer (both local)
  double answer = 0.0;
  double danswer = 0.0;

  // Constructing the equation
  for (int m = 1; m <=NInput; m++)
  {
    answer += ((15*Input[m])-5)*((15*Input[m])-5);
  }

  // Makes danswer equal to itself multiplied
  for(int n=1; n<=NInput; n++)
  {
    danswer += ((0.5*n)*((15*Input[n])-5));
  }

  answer += (danswer*danswer);
  answer += (danswer*danswer*danswer*danswer);

  Result.Input(answer, Input);

  nCall++;

  return (answer); 
}

double ZakharovFunction::NormValue(vector<double> &Input)
{
  // first check that point is not beyong any boundaries
  checkBoundary(Input);

  // Initialise vars answer and danswer (both local)
  double answer = 0.0;
  double danswer = 0.0;

  // Constructing the equation
  for (int m = 0; m <NInput; m++)
  {
    answer += ((15*Input[m])-5)*((15*Input[m])-5);
  }

  // Makes danswer equal to itself multiplied
  for(int n=0; n<NInput; n++)
  {
    danswer += ((0.5*(n+1))*((15*Input[n])-5));
  }

  answer += (danswer*danswer);
  answer += (danswer*danswer*danswer*danswer);

  Result.Input(answer, Input);

  nCall++;

  return (answer); 
}

// Calculates/builds the derivative
void ZakharovFunction::BuildDerivative(double *Input, double *Derivative)
{
  // NInput = dimension of function
  double mDeriv = 0.0;

  // Constructing the equation
  for(int j=1; j<=NInput; j++)
  {
    Derivative[j] = 0.0;
    mDeriv += ((0.5*j)*Input[j]);
  }

  for(int n=1; n<=NInput; n++)
  {
    Derivative[n] = (2*Input[n]) + (n*mDeriv) + ((2*n)*(mDeriv*mDeriv*mDeriv));
  }

  for(int mw = 0; mw<=NInput; mw++)
  {
    mDerivative[mw] = Derivative[mw];
  }
}

// used to test the function is working
void ZakharovFunction::test_function()
{
  ofstream out("output/test_of_zakharov.txt");
  
  out << "NUMBER_OF_DIMENSIONS: " << NInput << endl 
    << "      Subscript[x,k]      |     Value    |         Derivative[k]" << endl << endl;

  int NRun = 3;
  double *mInput;
  double *mDeriv;
  newA(mInput, NInput);
  newA(mDeriv, NInput);

  for(int run=1; run<=NRun; run++)
  {
    for(int n=1; n<=NInput; n++)
    {
      mInput[n]=(15.0*ran1()) - 5.0;
    }

    for(int k=1; k<=NInput; k++)
    {
      out << setw(11) << mInput[k];
    }
      // d[1] == 0, doesn't give derivative = 0
      // as d[1] is never 0.
      out << setw(16) << Value(mInput);

      BuildDerivative(mInput, mDeriv);

    for(int l=1; l<=NInput; l++)
    {
      out << setw(16) << mDerivative[l];
    }

    out << "\n\n";
    
  }

  deleteA(mInput);
  deleteA(mDeriv);
  out.close();
}

// used to print the solution (this is the old version, the new version is below)
void ZakharovFunction::print_sol(double *solution, double funcVal, char filename[])
{
  ofstream out(filename);	

  out << "Number of dimensions: " << NInput << endl;
  out << "Final Zakharov value = " << funcVal << endl;
  out << "Solutions in array I:" << endl;

  for (int i = 1; i <= NInput; i++)
  {
    out << solution[i] << "  ";
  }
 
  out << endl;

  out.close();
}

// used to print the solution to a file
void ZakharovFunction::print_sol(char filename[])
{
  ofstream out(filename);	
  
  // getBestValue() returns the best value stored in the tracker
  out << "Number of dimensions: " << NInput << endl;
  out << "Final Zakharov value = " << Result.getBestValue() << endl;
  out << "Solution coordinates:" << endl;
  
  // prints out the solution coordinates to a file
  Result.printSolution(out);

  out.close();
}

// distance from know global minimum which is located at the origin
double ZakharovFunction::distanceFromGlobalMinimum(double *solution)
{
  double toReturn = 0.0;

  for (int i = 1; i <= NInput; i++)
  {
    toReturn += solution[i]*solution[i];
  }

  double toReturnS = sqrt(toReturn);
  return toReturnS;
}